/****
This function returns a random number distributed as a power law
p(y)~y^gamma in the range ymin to ymax 


Written by MF&RD, Apr 2010.
****/


double ranpowerlaw(double ymin, double ymax, double gamma, gsl_rng * rand) {
  
  double x=gsl_rng_uniform (rand);
  return(pow(x*pow(ymax,gamma+1)+(1-x)*pow(ymin,gamma+1), 1.0/(gamma+1)));

}
